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Abstract 

We analyze two alternative methods for determining the exponent z of the contact process (CP) 
and Domany-Kinzel (DK) cellular automaton in Monte Carlo Simulations. One method employs 
mixed initial conditions, as proposed for magnetic models [Phys. Lett. A. 298, 325 (2002)]; the 
other is based on the growth of the moment ratio m(t) = (p 2 (t)) / (p(t)) 2 starting with all sites 
occupied. The methods provide reliable estimates for z using the short time dynamics of the 
process. Estimates of u\\ are obtained using a method suggested by Grassberger. 
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I. INTRODUTION 



The dynamics of spin models quenched from high temperature to the critical point is a 
subject of considerable current interest, because the initial phase of the relaxation (the so- 
called short time dynamics) carries important information about static and dynamic critical 
behavior. Using a phenomenological renormalization group analysis, Janssen, Schaub and 
Schmittmann |4j demonstrated the scaling law 

(M) (t) ~ m o t (1) 

characterized by the new critical exponent 9. (Here (•) denotes an average over initial con- 
figurations consistent with initial magnetization m , and over the noise in the stochastic 
dynamics.) Relation ([TJ holds for times t < t max ~ m z ^ x ° , where z is the dynamic ex- 
ponent, while mo denotes the initial magnetization, and Xq its anomalous dimension. The 
"critical initial slip" described by Eq. ((TJ) emerges from an initially disordered state, and is 
characteristic of a stochastic, far from equilibrium relaxation process. This approach also 
offers a way to determine the static critical exponents, since the scaling law used to deduce 
Eq. can also be applied to higher moments of the order parameter, yielding: 

(M k ) (t, mo) = (M h ) (L~ z t, L x °m ) (2) 

Models with absorbing states exhibit scaling behavior at the critical point marking the 
transition between active and absorbing stationary states, even though the stationary state 
is not given by a Boltzmann distribution The order parameter is the activity density 

p > 0, given by pit) = L~ d Yli=i a i ^ 0, where <jj = 1 (0) corresponds to the presence 
(absence) of activity at site % id is the dimensionality of the system). Activity is commonly 
associated with the presence of a "particle". For such models the following scaling law |l| 
has been conjectured: 

p(t)~r^f [{p-p c )t l ^^ z /N,p^ v » +e ] , (3) 

where po is the initial density, p is a temperaturelike control parameter, p c is its critical 
value, and N = L d is the number of sites. The exponent (3 is associated with dependence 
of the stationary value of p on the control parameter: p ~ (p — p^ 13 , while v\\ and z/j_ are 
the critical exponents associated with the correlation time (£n ~ \p — p c \ I ) and correlation 
length (£_|_ ~ \p — p c \~ u± ). Given the defining relation £y ~ £*, the dynamic exponent 
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z = u,./u±. Letting p — > p c and N — > oo, we expect the scaling function / in Eq. Q to 
have the property: 

u, u ~ 



f[0,0,u] 



(4) 

C, u — > oo 



where C is a constant. 

In a realization of the process at criticality (p = p c ) beginning with a completely filled 
lattice (p = 1), the activity density decays via the power law p(t) ~ t P,V W for t < 
L z . On the other hand, in realizations starting with only a single particle or active site 
(spreading process), the average number of particles increases with time: p(t) ~ tf 1 , where 
ri = (dv ± -2/3)/v n . 

An interesting crossover phenomenon in the evolution of p(t), between the initial increase 
(~ f) asymptotic decay (~ t ^"ll), emerges in a critical spreading process with low initial 
particle density. The crossover time t c is related to the initial density p y i a 

tc~Po ■ (5) 

A process starting with a single particle (p = 1/L) corresponds to p — > for large lattices, 
so that t c diverges and the spreading regime p(t) ~ t v extends over the entire evolution. 

A method for determining the critical exponent z, proposed in the context of the short 
time behavior of magnetic models |5j, suggests that we consider the function 

F 2 (t) = (6) 

d v- 

which has the asymptotic behavior i 7 ^) ~ i "n = t < z . Thus we can obtain z directly by 
analyzing short time simulation results with mixed initial conditions. 

In this paper we employ this approach to calculate z for the contact process (CP) and 
the Domany-Kinzel cellular automaton (DK), both known to belong to the universality 
class of directed percolation 2\. We should note that in the literature on absorbing-state 
phase transitions z is commonly used to denote the exponent governing the growth of the 
mean-square distance of particles from the original seed, in spreading simulations. To avoid 
confusion we denote the latter exponent as Z, so that R 2 ~ t z . The dynamic exponent is 
then related to the spreading exponent via z = 2/Z. 
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A method to estimate the exponent v\\ was suggested by Grassberger [ll|, who used the 
relation 

t l/uz ^ 



D(t)- dlnp 



d P P = Pc 

where in a simulation the derivative is evaluated numerically via 



fl(t ) 144^1, ( 8) 



2h \p(p c -h) 

which evidently requires data for values of p slightly off critical. In this context a reweighting 
scheme thatpermits one to study various values of p in the same simulation is particularly 
convenient [3j]. 

In the following section we present details on the models and our simulation technique. 
In section III we report and analyze the simulation results, and in section IV present our 
conclusions. 

II. MODELS AND SIMULATION METHOD 

The contact process was introduced by Harris as a toy model of epidemic propagation. It 
evolves in continuous time. Denoting the number of occupied nearest neighbors of site i by 
n i — a ji where (i) denotes the set of nearest neighbors of site i, the transition rates 

are 

p(0 — > l,n) — || (creation) 

(9) 

p(l — > 0, n) — 1 (annihilation). 
The model suffers a continuous transition at a critical value A c ; in one dimension A c = 
3.29785(8) U. 

nn 

As described in Refs. |2J, |3[ , our simulations employ a list of occupied sites and sample 
reweighting to improve efficiency Since the choice of sites in the dynamics is restricted to 
the occupied set, the time increment At associated with each event (annihilation or creation) 
is 1/Np, where N p is the number of occupied sites immediately prior to the event. A given 
realization of the process ends at a predetermined maximum time, or when all particles have 
been annihilated. Reweighting is used to study the effects of window size h in Eq. (jHj), in 
determining mi for the CP, using po — 1/L. 

The Domany Kinzel (DK) cellular automaton is a discrete-time process exhibiting a 
phase transition between an active and an absorbing phase of the same kind as in the 
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CP j8j. Each site of the lattice can be in one of two states, Oi = 1 (active) or a. t = 
(inactive). The transition probabilities for (Xj(t) given the values of its neighbors u±\{t — 1) 
are: P(1|0,0) = 0, P(1|1,0) = P(1|0,1) = p and P(l|l,l) = q. This model has a line 
of continuous phase transitions separating the active and absorbing phases in the p — q 
plane. For q = p(2 —p) the DK model is equivalent to bond directed percolation (DP), with 
p c = 0.644700. The critical behavior along the transition line falls in the DP universality 
class, with the exception of the point p = 1/2, q = I, corresponding to so-called compact 
DP 0,0. 



III. RESULTS 



We study the time-dependent order-parameter moments (p k (t)) in Monte Carlo simula- 
tions of the one-dimensional DK and CP, obtaining time series (pit)) po=1 / L and (p(£)) po=1 , 
for the two initial conditions described above. These quantities are combined to yield F 2 
defined in Eq.©. Initially we studied systems of L = 2048 sites in N s = 50000 independent 
realizations, each extending to t max = 1000. Analyzing our results, we find z = 1.5801(4) for 
the CP and z = 1.5804(2) for the DK. (For the CP data for times in the interval [30, 1000] 
are used; for the DK model the interval is [40,1000]. The exponent values are obtained 
from linear fits to the data on log scales.) Similar results are obtained for L = 4096. The 
evolution of F 2 for the two models is shown in Fig. 1. These results are in agreement with 
previous results on DP ( z = 1.580745(10) ^(|) obtained via a low-density expansion, and 
on the CP (z = 1.58077(2)) from exact diagonalization of the master equation jl^ . 
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FIG. 1: Time evolution of the F 2 to CP and DK. 



5 



To obtain our best estimate, we divide the time interval into subintervals equally spaced 
in Int. This avoids placing an undluy large weight on longer times, as would occur if each 
integer time were treated separate data point. 

In this case was used N s = 10000 samples for each independent seed, to a total of 5 
seeds, that gives the same N s = 50000 samples, however with lesser errors due the statistical 
correlations between samples. 

We study various intervals of lnt to obtain our estimates; for example, for [2.197, 6.904] 
we obtain 1/z = 0.62867(52) and 1/z = 0.62802(44) respectively for the CP and DK. 
Our best estimates are found respectively using the intervals [3.555, 6.904] for the CP and 
[4.9397, 6.904] for the DK, yielding 1/z = 0.63240(25) (z = 1.5813(6)) and 1/z = 0.63252(37) 
z = 1.5810(9)). These results are more realiable and are consistent with those of id] and 
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FIG. 2: Effective exponent z(t) versus 1/t to the CP and DK. Note the convergence to value 
z a; 1.58 in both cases. 

Numerically z t can be estimated from a least-squares linear fit; it is plotted versus t" 1 , t a 
being the geometric mean of the t— values over the subintervals. The exponents z t and vn 
are similarly obtained, and an extrapolation t^ 1 — > is performed. In Fig 2 we show this 
extrapolation for the CP and DK, showing convergence towards z ~ 1.58, that reflects the 
universality between the two models. 

The ratio 

m(t ) = <f«$f=L (10 ) 
converges to the expected value for models in the DP universality class, ~ 1.1174 to 



DP in (1 + 1) dimensions [13]. At short times, m — 1 = var(p)/ (p) z increases as a power 
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law. The associated exponent is found by noting that 



X = L d var(p) ~ t*g[(p - p c pt} ^ (p - ptf (11) 

so that g(x) ~ a; -7 ^ for large x. Observing that (j) = y/f\\ = (dv± — 2/3) /j/n, and that 
(p) ~ £ II , we expect m — 1 ~ 

t*A. For DP in (1 + 1) dimensions, 1/z = 0.6326. Our 
simulations of the CP (L = 5000, N s = 10000), give z~ x = 0.634(3). 




FIG. 3: Plot of time evolution of D{t) in the CP and DK. 

In the determination of v\\ via Eq. ((7j), we used t max = 2980. In Fig. 3 we show the 
derivative D(t), and in Fig. 4 the effective exponent v\\{t) is plotted versus 1/t for the 
CP and DK, for the case p — 1- We find v\\(t) ~ 1.70, as t — » oo. We refined this 
estimate using lattices of L = 2048, 4096 and 8192 for the DK model, to gauge finite-size 
effects. The estimates for u\\ (using a total of N s = 10000 realizations, and Ah = 25 = 
0.002) were 1.623(2), 1.625(3) and 1.648(5) respectively. These differ from the accepted 



value v\\ = 1.733847(6) [2j, |10|. Reducing Ah by an order of magnitude, to 0.0002, we 
obtain uu = 1.731(9), consistent with the expected result. In the CP simulations we used 
the reweighting method proposed in j3j. The sample of realizations generated at A c is 
reweighted for nearby values, A c ± n8. Using 5 = 0.001, we determine v\\(t) (figure 5) to 
study the influence of Ah. The last ten points of each curve are extrapolated to 1/t — *■ 
to estimate v\\\ as 5 is reduced, our estimate approaches the expected value uu ~ 1.73. 
(For 5 < 0.001 the curves become indistinguishable.) The value found via extrapolation is 
v\\ = 1.734(8), in agreement with the accepted value. 
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FIG. 4: Effective exponent versus 1/t to the CP and DK with po = 1. 




0.001 0.002 0.003 0.004 0.005 0.006 0.007 

1/t 



FIG. 5: Effects of Ah size in the curves of effective exponent in the CP modelo to po = 1/L. 
Non-equilibrium reweighting and the continuous time MC method were used. 

A interesting way of measure the direct influence of Ah in u\\ can be seen in a plot of u\\ 
as function of Ah. We note the expected behavior of v\\ — > 1.73.... (figure 6). 

IV. SUMMARY AND CONCLUSIONS 



We apply several simulation methods based on analysis of short time simulations to deter- 
mine the critical exponents z and v» in models exhibiting a continuous phase transition to an 
absorbing state. The ratios F 2 (t) = (p) po=1/L d/ (p) 2 po=l and m(t) = (p(t) 2 ) po=1 / {p{t)f po=1 
have shown to be useful in alternative methods to determine critical exponents of those mod- 
els. Our estimates for the dynamic critical exponent z (1.5813(6) for the CP and 1.5810(9) 
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FIG. 6: Plot of z/ii as a function of n, where Ah = 2nd. 

for the DK cellular automaton) are in good agreement with accepted results in the literature 
obtained by other techniques such as low-density expansion and exact diagonalization of the 
master equation. 

In order to improve the efficiency in determining the critical exponent v\ \ we employ a list 
of occupied sites and sample reweighting. We also check the influence of h, the increment 
used in evaluating the numerical derivative in Eq. (jHJ). The estimates for uu appear to be 
very sensitive to increment size. Using Ah = 25 = 0.0002 we obtain our best estimates 
which are v\\ = 1.731(9) for the DK cellular automata and v\\ = 1.734(8) for the CP process. 
These results are in agreement with, though considerably less precise than, the best estimate 
in the literature, u\\ = 1.733847(6) [lOj. We believe that the methods investigated here will 
be useful in the analysis of other models with absorbing states, in particular, in establishing 
the universality class using short time simulations, which are typically less computationally 
demanding than studies of the stationary process. 
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